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Recently, Li et al. (Bioinformatics 27(19), 2686-91, 2011) proposed a method, 
called Differential Equation-based Local Dynamic Bayesian Network (DELDBN), 
for reverse engineering gene regulatory networks from time-course data. We 
commend the authors for an interesting paper that draws attention to the 
close relationship between dynamic Bayesian networks (DBNs) and differen- 
tial equations (DEs). Their central claim is that modifying a DBN to model 
Euler approximations to the gradient rather than expression levels themselves 
is beneficial for network inference. The empirical evidence provided is based 
on time-course data with equally-spaced observations. However, as we discuss 
below, in the particular case of equally-spaced observations, Euler approxima- 
tions and conventional DBNs lead to equivalent statistical models that, absent 
QO artefacts due to the estimation procedure, yield networks with identical inter- 

(T) gene edge sets. Here, we discuss further the relationship between DEs and 

C*~) conventional DBNs and present new empirical results on unequally spaced data 

— which demonstrate that modelling Euler approximations in a DBN can lead to 

improved network reconstruction. 

A dynamic Bayesian network (DBN) is a Bayesian network (i.e. a graphical 
' . I model based on a directed acyclic graph) with an explicit time index. Consider 

^ biochemical time-series expression data Xi(t) where i € {1 . . .p} indexes genes 

(or other molecular variables of interest) and t € {1 . . . T} indexes time. In 
the DBNs considered in Li et al. ( 2011[ ) (with linear Gaussian conditionals and 
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03 edges only from one time slice to the next) mean expression at time t + 1 is 

modelled as a linear function of expression at time t 

X^t + V-N^^X^t),*!] (1) 

where /3y is a parameter describing the influence of gene j on target i, Oi is 
a noise parameter and N(fj,, v) denotes a Gaussian distribution with mean \i 
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and variance v. For convenience we refer to models of the form in Eqn. [T] as 
"conventional" DBNs. 

In contrast, consider a DE that models the gradient as a deterministic linear 
function with parameters fy: 



dXj{t) 
dt 



p 



DELDBN combines DBNs (Eqn. [T| with DEs (Eqn. [2} by using an Euler 
gradient approximation 



Xjjt + l) -Xj(t) 
At 



N Y>i*i(*).°? ( 3 ) 



where Af denotes the interval, in units of time, between observations with time 
indices t and t + 1 . 

In the terminology of regression, the left hand sides of Eqns. [T] and [3] are 
responses that are modelled using predictors Xj , with independent samples in- 
dexed by t. The main claim of Li et al. is that improved performance in network 
reconstruction may be achieved by modelling the response as the Euler gradient 
(Eqn. [3]) rather than the observed value (Eqn. [TJ of gene expression at a given 
time, provided the time interval At between samples is not too large. Eviden- 
tial support for this conclusion is provided using data from a synthetic gene 
regulatory network "IRMA" that was constructed in Saccharomyces cerevisiae 



( |Cantone et q/.[|2009[ ) 



DELDBN carries out inference regarding network topology using Markov 
blankets, facilitated by a heuristic search implemented in the R package BNLearn 



(Scutari 2010). The Markov blanket for a node i in a Bayesian network is the 
set of nodes comprising i's parents, its children, and its children's other parents. 
For the DBNs considered in Li et al. the Bayesian network is bipartite and 
a target i has no children, only parents (as depicted in Fig. 1 in Li et al.). 
Therefore the Markov blanket MB(i) is identified with the parents of i in the 
DBN. According to the Euler gradient model (Eqn. |3| the Markov blanket of a 
target Xi is given by the non-zero coefficients 

MB Eu]er (Xi) = {Xj : 0}. (4) 

When the time At between samples is constant, the Euler gradient model is 
simply a reparameterisation of the conventional DBN: 

p p 

Xi(t + 1) = Xi (t) + AtJ2 hjXj (t) = E ~^ X i W ( 5 ) 
j'=i j=i 

AtPij 1 £ j 
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In this case of equally spaced observations, under the conventional DBN the 
Markov blanket is 



MBc-dbnPQ) = {Xj : hj + 0} (7) 

which differs from Eqn. [3] only in the possible presence/absence of a self-loop 
Xi, since for i ^ j we have that ^ if and only if /Jy ^ 0. Thus, for 
equally spaced observations the Markov blanket is invariant (up to inclusion 
of a self-loop) to the reparameterisation that occurs in going from a conven- 
tional DBN to Eulcr gradient responses. This means that the two formulations 
are equivalent with respect to inference regarding the inter-gene edge set when 
equal sampling intervals are used. Therefore in order to distinguish between 
the reverse-engineering performance of these approaches, it is essential to con- 
sider the regime in which data are sampled unevenly in time. However, results 
reported by Li et al. were obtained using only data sampled evenly in timeQ 

Nevertheless, we fully agree that the relationship between dynamics and 
DBNs merits careful investigation. In particular many time-course datascts in 
bioinformatics are obtained with unequal sampling intervals. Then, the equiv- 
alence between conventional DBNs and Euler gradient models does not hold, 
making the choice of formulation an important question. We therefore under- 
took empirical comparison of a range of modelling approaches, which may all be 
viewed as variations of the well studied variable-selection problem in linear re- 



gression (Oates and Mukherjee 2011). This subsumes both conventional DBNs 
and the Euler gradient model discussed above. 

In order to obtain unevenly sampled data from the IRMA network studied 



in Li et al., we used the differential equation (DE) model described by Cantone 



et al. (20091. This model has been demonstrated to provide good fit to the 
IRMA data. We generated data at unevenly spaced times 0, 1, 2, 4, 6, 10, 15, 
20, 25, 30, 40, 50, 60, 80, 100, 140, 180, 220 and 280 minutes, adding Gaussian 
measurement error with variance set to give a signal-to-noise ratio equal to 20. 
A typical dataset generated in this way is presented in Figure [T] 

We carried out inference within a Bayesian framework. Denote by y the 
responses (for target i, for simplicity we suppress dependence on i in what 
follows), either as in a conventional DBN or Euler approximations. Let D 
denote the design matrix constructed according to the Markov blanket MB (for 
notational simplicity we leave dependence on MB, i.e. the graph structure, 
implicit below), and let m be the number of columns in D. Assuming additive 



1 lt is interesting to ask why the output of DELDBN on the IRMA data appeared to 
improve using Euler gradient responses. The growth-shrink (GS) algorithm (Margaritis and 
|Thrun[ |1999[ l was used to infer Markov blankets from data. GS proceeds by carrying out 
conditional independence tests, based in this case on the Pearson correlation coefficient r 
and statistic rw ^Z^z ~ Tt—2, where 7t-2 denotes the t-distribution with T — 2 degrees 
of freedom. This particular approximate approach to Markov blanket identification while 
computationally efficient is not invariant to the reparametrisation relating the conventional 
DBN and Euler gradient models. Since these models are structurally equivalent with respect 
to inter-gene edges, this suggests that the improved performance of DELDBN on inter-gene 
edges reported in Li et al. may be an artefact due to the specific estimator used. 
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Figure 1: Simulated data, typical dataset. The DE model of |Cantone et at] 
(2009) was used to simulate data from the IRMA network at uneven time in- 
tervals. 



Gaussian error we get y = Dj3 + e, e ~ N(0, <J 2 I) where / denotes the identity 
matrix, f3 collects together all coefficients and a 2 is the variance of the error. 
We score models by integrating the corresponding likelihood against a prior for 
(/3,cr 2 ). Here, we used a g-prior ( |Zellner[ [l986| /3 | a 2 ~ N(0,a 2 T(D'D)- 1 ) 
for coefficients and 7r(cr 2 ) cx 1/cr 2 , where n is the sample size in the regression 
sense. This leads to closed-form marginal likelihood 



w{y | MB) 



f 



1 + T 



m/2 r 



y v 



T 



1 + T 



n -(T-l)/2 



y y 



(8) 



where y — D(D'D) 1 D'y. Network inference is carried out by Bayesian model 
averaging, using the posterior probability 

{Xj € MB(i)}7r(y | MB)tt(MB) 



P(j regulates i) = 



MB 



Emb- <y I mb')tt(mb') 



(9) 



to score a directed edge from gene j to target i, where I{^4} = 1 is A is true, 
otherwise I{A} = 0. 

In experiments below we take a network prior which, for each target i, is 
uniform over the number of predictors m, up to a maximum permissible in- 
degree d max , that is 7r(MB) cx \\ i ( ^ ) I{m, < d max }, but note that richer 
network priors are available in the literature (Mukherjee and Speed |2008 ). A 
Markov blanket estimator is obtained by thresholding posterior edge probabil- 
ities; for threshold r this gives a network with estimated edge-set E = {(j, i) : 
P(j regulates i) > r}. For small maximum in-degree <i max , exact inference by 
enumeration of variable subsets may be possible. Otherwise, Markov chain 
Monte Carlo (MCMC) methods can be used to estimate posterior edge proba- 
bilities ( Ellis and Wong| |2008 Friedman and Koller, 2003). In the experiments 
here we use exact inference by enumeration, with <i max = 2. 
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Figure 2: IRMA network, area under receiver operating characteristic curves 
(AURs). AURs were calculated for a range of modeling approaches, based on 
data simulated from the IRMA network at uneven time intervals. Key: "Resp" 

- [Euler] Euler derivative approximations / [C-DBN] conventional DBN; "Pred" 

- [Prod] using products XiXj of predictors / [Std] otherwise; "Lag" - [Yes] 
additional lagged predictors / [No] otherwise. 



Under the framework outlined above, we assessed network reconstruction 
using both conventional DBNs and Euler gradients. Since the DE model of 



Cantone et al. (2009) is nonlinear, it is natural to also investigate whether the 



use of products XjXk of predictors (in addition to linear predictors) improves 
network reconstruction, by capturing some nontrivial aspects of the dynamics. 
Similarly, as suggested in the discussion of Li et al, since the DE model of |Can-| 



tone et al. ( 2009 1 contains a delay term, it is interesting to investigate whether 



the use of lagged predictor variables improves performance. The approaches we 
considered can be summarised as: 



Predictor set 
Lagged predictors 
Response 



{ Standard, Product } 
{ No, Yes (lag « T/10) } 
{ Conventional DBN, Euler gradient} 

Performance was assessed using the area under receiver operating characteristic 
curves (AUR), equivalent to the probability that a randomly selected true edge 
has a higher score than a randomly selected false edge; higher values of AUR 
correspond to better performance. We carried out inference for 1000 sampled 
datasets for each method to obtain distributions over AUR scores, as shown 
in Figure [2] Inference based on the Euler gradient response outperforms the 
conventional DBN, supporting the central claim of Li et al. The use of products 
of predictors together with the inclusion of lagged predictors led to slightly 
improved performance. 



5 



In summary, we presented empirical evidence that Euler approximations to 
dynamics coupled with DBNs can be useful in reverse engineering gene regula- 
tory networks. Furthermore we showed how such models may be viewed in a 
regression framework, for which there exists a wide literature on variable selec- 
tion. However our investigation was somewhat idealised, since in practice data 
are often obtained under destructive sampling and averaging over large numbers 
of cells. An extended discussion on the relationship between cellular dynamics, 
nontrivial observation processes and linear regression may be found in |Oates 



and Mukherjee (20111. 
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